%!TeX encoding = utf-8
%!TeX program = xelatex
\documentclass[12pt,a4paper]{article}                                                                            
\usepackage{amsmath,amsfonts,amsthm,amssymb}%AMS数学符号字体等的宏包
\usepackage{amsrefs}                           %ams参考文献宏包, 不需要的用户请删除此行!
%\usepackage{hyperref}                          %生成超链接
%\usepackage[colorlinks,linkcolor=black,anchorcolor=black,citecolor=black]{hyperref}

\usepackage{geometry}                        %页面设置
\usepackage{xeCJK}                               %中文字体宏包
\usepackage{graphicx}                          %插图
\usepackage{setspace}  
\usepackage{caption}
\captionsetup{font={small}}
\usepackage{listings} %C, C++等代码高亮，需要配合xcolor宏包支持
\usepackage{tabularx}
\usepackage{xcolor}

%%%%%%%%%%实用功能%%%%%%%%%%%%%%%%%%%%%%
\usepackage{indentfirst} % 解决 CJK 手段缩进问题,\setlength{\parindent}{2em}



\usepackage{cite}

%%%%%%%%%页面设置%%%%%%%%%%%
\setlength{\textwidth}{162mm}
\setlength{\textheight}{231mm}
\setlength{\headheight}{0cm}
\setlength{\topmargin}{-0.1cm}
\setlength{\oddsidemargin}{0cm}
\setlength{\evensidemargin}{0cm}
\setlength{\parskip}{1mm}
\setlength{\unitlength}{1mm}
\setlength{\parindent}{2.06em}
\linespread{1.5}

%%%%%%%%%%重命名%%%%%%%%%%
\renewcommand\refname{参考文献}
\renewcommand\contentsname{目录}
\renewcommand{\figurename}{图}
\renewcommand{\tablename}{表}
\renewcommand\abstractname{\bf\large 摘~~要}


%%%%%%%%%%定义编号规则%%%%%%%%%%
\renewcommand{\thefigure}{\arabic{section}.\arabic{figure}}
\renewcommand{\thetable}{\arabic{section}.\arabic{table}}
\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}


%%%%%%%%%%CJK字体设置%%%%%%%%%%
\setCJKmainfont[BoldFont=STSongti-SC-Bold, ItalicFont=STKaiti]{STSong}
\setCJKsansfont[BoldFont=STHeiti]{STXihei}
\setCJKmonofont{STFangsong}

\newcommand\dd{\mathrm{d}}

\begin{document}

% Listing环境全局设置
\lstset{
%numbers=left,
%numberstyle=\tiny,
keywordstyle=\color{blue!70}, 
commentstyle=\color{red!50!green!50!blue!50},
frame=shadowbox,
rulesepcolor=\color{red!20!green!20!blue!20},
escapeinside=``,
xleftmargin=2em,
xrightmargin=2em,
aboveskip=1em
}





\def\chntoday{\the\year~年~\the\month~月~\the\day~日}

\title{并行求解三维热传导方程}
\author{姚伟鹏\footnote{电话：13370130802，邮箱：ywp@pku.edu.cn，学号：1501110118 }}
\date{\chntoday}

\maketitle                            % 生成标题

\tableofcontents                      % 插入目录

\newpage

\section{问题模型}

\subsection{数学描述}

考虑如下热传导问题：		
\begin{eqnarray}
	\left \{
	\begin{array}{lll}
		\displaystyle \partial u/\partial t = \Delta u \qquad\quad in \quad \Omega \\
		\displaystyle u|_{t=0} = 1/2 \qquad\quad in \quad \Omega \\
		\displaystyle u = \left \{
		\begin{array}{l}
			1 \qquad\qquad \vec{r} \in \quad \Gamma_{out} . \\
			0 \qquad\qquad \vec{r} \in \quad \Gamma_{in} . \\
		\end{array}
		\right.
	\end{array}
	\right.
\end{eqnarray}

其中，各个区域分别为：
\begin{eqnarray}
	\left \{
	\begin{array}{cr}
		\Omega = [-2,2]^3 \backslash \{|x| + |y| + |z| < 1\} \\
		\Gamma_{out} = \{|x| = |y| = |z| = 2\} \\
		\Gamma_{in} = \{|x| + |y| + |z| = 1\} 
	\end{array}
	\right.
\end{eqnarray}
整个区域如图\ref{Grid}所示。
\begin{figure}[h!]
	\center{
  	\includegraphics[width=\textwidth]{files/grid.pdf}}\\
  	\caption{求解区域}
	\label{Grid}
\end{figure}


根据题意，上述模型的定解问题为初边值问题：

（1）当 $t=0$时，指定初始条件为$u(x,y,z,0)=u_0(x,y,z)\equiv1/2$.

（2）对于区域$\Omega$的外围正方体边界$\Gamma_{out}$，指定边界条件：
	\begin{equation}
	\label{bc_out}
		u(x,y,z,t) = 1,(x,y,z) \in \Gamma_{out}, t \geq 0
	\end{equation}
	
（3）对于区域$\Omega$的内部正八面体边界$\Gamma_{in}$，指定边界条件：
	\begin{equation}
	\label{bc_in}
		u(x,y,z,t) = 0,(x,y,z) \in \Gamma_{in}, t \geq 0
	\end{equation}
	

\subsection{物理预期}

从最简单的角度去考虑上述模型所代表的物理意义，即：

一个初始温度均匀分布的立方体块，在$t>0$时，放置于具有极其高温的外围环境中，使其外边界温度瞬间变为归一化后的$1$。
然而，由于某种原因，在内部的正八面体表面温度骤降为0。

不难预测，随着时间的推移，高温区会从外部边界逐渐向内蔓延；而在内部的正八面体内，则是低温区会从正八面体表面逐渐向内蔓延。
最终形成一个内部为零度且绝热的正八面体，被包裹在外部高温的正方体内的分布。

这是我们检验程序运行正确的标准之一。

\subsection{求解策略}

我们将采用有限差分法来求解热传导方程，首先对求解区域进行网格剖分，然后构造相应的数值计算格式，包括在$\Omega$区域内网格点的有限差分格式，
以及在$\Gamma_{in}$与$\Gamma_{out}$上的边界处理（这里由于两个边界均为简单的Dirichlet边界，故只需直接对边界点赋值即可）。

具体求解过程见第二章。

\section{数值求解策略}

\subsection{网格离散}
由于区域$\Omega$具有对称性，可以采用结构网格剖分，直接选用正方体作为网格单元，剖分节点$P_{i,j,k}$为：
\[ P_{i,j,k} = (ih,jh,kh),|i| +|j| +|k| \leq  2N\]
其中，$h$为空间步长，$Nh=1$。

与节点$P_{i,j,k}$相连接的六个节点（假设也在$\Omega$内）为：
$P_{i-1,j,k}$, $P_{i+1,j,k}$, $P_{i,j-1,k}$, $P_{i,j+1,k}$, $P_{i,j,k-1}$, $P_{i,j,k+1}$。

节点在边界上的充要条件是：
\begin{eqnarray}
    \label{BC1}
    P_{i,j,k} \in \Gamma_{in} \Leftrightarrow |i| +|j| +|k| = N \\
    \label{BC2}
    P_{i,j,k} \in \Gamma_{out} \Leftrightarrow |i| = |j| = |k| = 2N
\end{eqnarray}

\subsection{热传导方程的差分格式}

在区域$\Omega$内，利用有限差分法来求解热传导方程。
记$u_{i,j,k}^n = u(ih,jh,kh,n\tau)$，其中$\tau$为时间步长。

引入微分算子：$D_x = \partial / \partial x$，$D_y = \partial / \partial y$，$D_z = \partial / \partial z$，
则$\Delta = D_x^2 + D_y^2 + D_z^2$。
利用差商代替微商的方法，即用适合的一阶差商的组合来替换方程中出现的各阶微商，从而构造相应的差分格式。

构造差分格式的工具之一是差分算子，其中最基本的是一阶中心差分算子，如下所示：
\[ \delta_x u_{i,j,k} = u_{i+\frac{1}{2},j,k} - u_{i-\frac{1}{2},j,k}\]
连续应用两次一阶中心差分算子$\delta_x$，得到常用的二阶中心差分算子：
\[ \delta_x^2 u_{i,j,k} = \delta_x(\delta_x u_{i,j,k}) = u_{i-1,j,k} - 2u_{i,j,k} + u_{i+1,j,k}\]
同理可得：
\[ \delta_y^2 u_{i,j,k} = \delta_y(\delta_y u_{i,j,k}) = u_{i,j-1,k} - 2u_{i,j,k} + u_{i,j+1,k}\]
\[ \delta_z^2 u_{i,j,k} = \delta_z(\delta_z u_{i,j,k}) = u_{i,j,k-1} - 2u_{i,j,k} + u_{i,j,k+1}\]
用差分近似可得：
\begin{equation}
	\label{fdm}
   	 \frac{u_{i,j,k}^{n+1}-u_{i,j,k}^n }{\tau}
    	= \frac{ \delta_x^2 u_{i,j,k}^n+ \delta_y^2 u_{i,j,k}^n+  \delta_z^2 u_{i,j,k}^n }{h^2}
\end{equation}
容易看出这个格式的截断误差为$O(\tau + h^2)$,而且当$\tau \leq h^2/6$时是数值稳定的。

\subsection{处理初值与边界条件}

上面的数值格式只在区域的内点有效，即所用到的数据必须是来自于内点的数据。
从网格节点的分布可以看出，计算$u^{n+1}$的内点的数据，只会用到$u^n$的内点和边界点的数据。
对于$t=0$的数据，可以直接从初值条件得到，即：
\begin{equation}
	\label{initial}
	u_{i,j,k}^0 = 1/2
\end{equation}
在边界点处，需要单独处理边界处的计算格式。

由于在我们的问题中，两个边界条件均为简单的Dirichlet边界，故只需直接对边界点赋值即可。

（1）当$P_{i,j,k} \in \Gamma_{in}$时，可以直接赋值：
\begin{equation}
	\label{Boundary1}
	u_{i,j,k}^n = 0
\end{equation}

（2）当$P_{i,j,k} \in \Gamma_{out}$时，可以直接赋值：
\begin{equation}
	\label{Boundary2}
	u_{i,j,k}^n = 1
\end{equation}


%\section{差分格式}

\section{MPI并行设计}
MPI并行算法必须考虑两个关键问题, 即数据分块和数据通信, 合理的数据分块和数据通信将直接影响到MPI并行的性能。

\subsection{分组原则}
事实上，在我们这里的问题中，由于计算$u_{i,j,k}^{n+1}$仅仅能用到$u_{i,j,k}^n$周围少量节点处的数据，
所以可以直接将节点做分组，每一组用一个进程来做计算，即可实现并行。

\subsection{并行模式}
假设有$M$个进程，若采用主从模式，一个进程充当根进程，负责处理数据的分发、收集以及读写；其他的进程则负责计算。
于是可以将节点分为$M-1$组。虽然每组的节点数不一定相同，但原则上尽量保证每组的节点数平均分配，从而保证进程的负载平衡。

\subsection{数据通讯}
实现数据的传输有两种办法：聚合通信与单边通信。

聚合通信需要根进程先将上一时刻的数据分发给其他进程，每一个进程完成自己的计算任务之后，再由根进程进行收集，这样即实现了一个时间步的计算。
与之不同，单边通信则需要给每一个进程分配窗口，在计算过程中根据需要读取和写入数据。

显而易见，由根进程收集和分发数据的传输模式会给根进程带来较大的压力。
为了分担根进程收发数据的压力，可以在子进程之间传输数据，以实现数据的共享，即使用单边通信的方式。

（1）在数据分发阶段，最简单的方法是将上一时刻的所有数据广播到每一个进程，这个方法实现比较简单，但效率比较低。
因为每个进程在计算时并不需要所有的数据。

一个改进的办法是只分发那些计算时需要的数据，这首先要求我们识别出哪些数据是需要的，而且还要想办法使得这些数据的调用较为方便的实现。
为了实现上述操作，如图\ref{Group}所示，可以建立一个子数据数组。
它是原始数组中的某一段，虽然计算时并不一定会用到这一段中的所有数据，但使用该子数组仍旧比发送整个数组要高效一些。
不难发现，这个过程还可以通过相连接的几个进程之间点对点通讯做数据交换来实现。
\begin{figure}[h!]
	\center{
  	\includegraphics[width=\textwidth]{files/Group.jpg}}\\
  	\caption{每个进程上用到的数据}
	\label{Group}
\end{figure}

（2）在数据的收集和输出阶段，相对比较容易实现。
因为收集的数据刚好是连续且不重叠的，而且只需要在根进程上做输出即可。


\section{数值计算}


\subsection{数据结构}
在第二章中，我们给出的数值计算格式虽然是明确地，但是在计算过程中仍需要明确数据的储存和调用方式。
数据主要包括网格信息和温度计算结果，这里都蚕蛹数组形式储存。

网格信息包括节点的编号，坐标，以及节点之间的连接信息。

（1）首先考虑节点的编号。
每一个节点都有一个唯一的编号$(i,j,k)$，它和节点的坐标可以直接对应起来。
然而，在储存和调用时最简单的方式是单指标的形式，所以这里我们要建立一个节点编号表，从而将原来的三指标编号转换为单指标编号。
作为一个二维数组，该节点编号表的结构如表\ref{PointsIndex}所示。
其中，最后一列用来标记节点是属于哪一个区域或哪一个边界。

\begin{table}[h!]
	\centering
	\caption{节点编号表(Points Index)}
	\label{PointsIndex}
	\begin{tabular}{|c|c|c|c|c|}
		\hline
		节点编号 & $i$ & $j$ & $k$ & 节点位置\\
		\hline
		0 & $i_0$ & $j_0$ & $k_0$ & $\Omega$ \\
		\hline
		1 & $i_1$ & $j_1$ & $k_1$ & $\Gamma_{in}$\\
		\hline
		2 & $i_2$ & $j_2$ & $k_2$ & $\Gamma_{out}$\\
		\hline
		... & ... & ... & ....  & ...\\
		\hline
	\end{tabular}
\end{table}

（2）接下来考虑节点的连接信息。 
同样作为一个二维数组，节点连接信息的结构如表\ref{PointsConnectTable}所示
（取$\Omega$区域的最右前下角处的节点编号为$0$，如图\ref{pct}所示）。
其中，$-1$表示该节点不在计算区域内。
\begin{figure}[h!]
	\center{
  	\includegraphics[width=\textwidth]{files/pct.pdf}}\\
  	\caption{节点的连接信息}
	\label{pct}
\end{figure}

\begin{table}[h!]
	\centering
	\caption{节点连接表(Points Connection Table)}
	\label{PointsConnectTable}
	\begin{tabular}{|c|c|c|c|c|c|c|}
		\hline
		节点编号 & 左节点 & 右节点 & 前节点 & 后节点 & 上节点 & 下节点\\
		\hline
		0 & -1  & 1 & -1 & 2 & 3 & -1\\
		\hline
		1 & 0 & 4 & -1 & 5 & 6 & -1\\
		\hline
		2 & -1 & 5 & 0 & 7 & 8 & -1\\
		\hline
		3 & -1 & 6 & -1 & 8 & 9 & 0\\
		\hline
		... & ... & ... & ....  & ... &... &...\\
		\hline
	\end{tabular}
\end{table}

（3）最后，考虑温度的计算结果。
各个时刻各个节点处的温度值，本可以用一个较大的二维数组来储存。
但是由于计算下一时刻的温度仅需要用到当前时刻的数据，为了节省内存，
我们在这里只将之前的数据以文本的形式记录下来。
以此，只需要用两个一维数组来储存当前时刻和下一时刻的温度值，数组的长度均为节点总数。


\subsection{生成网格数据}
考虑在之前建立的网格剖分策略下，节点的编号表和连接表如何确定。

（1）首先分析整个计算区域的节点数目，因为$Nh=1$，而区域$\Omega$所在的立方体为$[-2,2]^3$，
故节点的数目为$A_N=(4N+1)^3$。

（2）接下来，对节点进行单指标编号。这里使用字典排序法实现：
\begin{lstlisting}[language={[ANSI]C++}]
int AN = 0;
for(int i=-2*N; i<=2*N; i++)
    for( int j=-2*N; j<=2*N ;j++)
	for( int k=2*N; k<=2*N; k++ )
	    // 第 AN 个节点为(i,j,k)
	    Table[AN++] = (i,j,k); 
\end{lstlisting}

（3）节点所在区域或边界可以使用公式（\ref{BC1}）和公式（\ref{BC2}）来确定。

（4）并行时节点的分组。要将所有$A_N$个节点分为$M-1$组，可将编号连续的若干个节点分为一组，
每组的节点数尽量接近，用一个一维数组来记录每一组的起始和结束编号。

\subsection{计算步骤}

我们采用显示格式来计算，算法步骤如下：

（1）确定网格。

给定剖分网格参数$N$之后，确定网格节点，并对节点进行编号，给出节点之间的连接信息，以及节点的所属区域或边界的信息。

（2）初始化。直接按照公式（\ref{initial}）赋值即可。

（3）根据进程数对网格节点进行分组。

（4）数据交换，使每个进程拿到需要的数据段。

（5）在每个节点处计算下一时刻的温度：

如果节点在$\Gamma_{in}$上，那么直接按照公式（\ref{Boundary1}）进行赋值；
如果节点在$\Gamma_{out}$上，那么直接按照公式（\ref{Boundary2}）进行赋值；
如果节点在$\Omega$的内点上，那么按照差分公式（\ref{fdm}）进行计算。

（6）数据收集到根进程，输出。

（7）结束程序。

这里，由物理预期不难得知，理论上当时间充分大之后，温度分布的变化趋于平缓，这时继续计算也没有太大意义。
所以我们设置计算终止的条件是：
\begin{equation}
	\label{stop}
	\parallel u(x,y,z,n\tau)-u(x,y,z,(n-1)\tau)\parallel_{L^\infty} \leq \varepsilon
\end{equation}
在随后的具体计算算例中，我们取$\varepsilon=10^{-9}$。


\section{结果与分析}

利用上面设计的算法，可以得到物体的温度分布变化。
为了验证计算结果的准确性，我们从计算精度和物理定性分析两个方面对结果做了分析。

本文使用的并行环境是MPICH-3.2。
CPU信息：Intel(R) Xeon(R) CPU E5-2650 @ 2.00 GHz $\times$ 16。

\subsection{计算结果}
这里我们选取$N=10,\tau=0.9\times \frac{1}{6N^2}$作为演示结果的算例参数。
此时，节点的数目为68921个，计算10000步之后，连续1000步的温度差小于$\varepsilon = 10^{-11}$，
此时的温度分布如图\ref{dist_3D}和\ref{dist_2D}所示。
\begin{figure}[h!]
	\center{
  	\includegraphics[width=\textwidth]{files/u_dist_3D.jpg}}
	%\includegraphics[width=0.5\textwidth]{files/dist_2D.pdf}}
  	\caption{温度三维空间分布（八分之一的$\Omega$区域，即\{$x>0 \cap y>0 \cap z>0$\}）}
	\label{dist_3D}
\end{figure}

\begin{figure}[h!]
	\center{
  	%\includegraphics[width=0.5\textwidth]{files/dist_3D.pdf}}
	\includegraphics[width=\textwidth]{files/u_dist_2D.jpg}}
  	\caption{温度二维空间分布（x-y平面）}
	\label{dist_2D}
\end{figure}

\subsection{物理定性分析}
从上面的温度分布可以看出，最高温度分布在$\Gamma_{out}$边界上，最低温度分布在$\Gamma_{in}$边界及其内部，
中间区域的温度由内向外逐渐升高。

这与我们在第一章中的物理预期基本一致。

%\subsection{计算精度验证}

\subsection{并行效率分析}
为了分析我们设计的并行算法的效率，这里给出了不同计算规模$N$，不同进程数$M$时，计算50000步所需要的时间比较，如表\ref{time_needed}所示。
\begin{table}[h!]
	\centering
	\caption{计算50000步所需的运行时间(s)}
	\label{time_needed}
	\begin{tabular}{|c|c|c|c|c|c|c|}
		\hline
		N & M=2 & M=4 & M=7 & M=10 & M=13 & M=16\\
		\hline
		10 &  25.9  & 10.9 & 6.6 & 5.4 & 4.0 & 2.9\\
		\hline
		20 & 216.7 & 98.9 & 59.4 & 42.4 & 32.1 & 25.4\\
		\hline
	\end{tabular}
\end{table}

据此，可以计算出并行加速比和并行效率：

\begin{center}
加速比 = 串行执行时间 / 并行执行时间

并行效率 = 串行执行时间 / (并行执行时间 $\times$ 进程数)
\end{center}

图\ref{accel}为取$N=10$与$N=20$时并行计算加速比和效率曲线图。
由于我们的程序是在每个节点拥有16个核的Cluster上运行的，所以当$M \leq 16$时，
一方面，加速比都是随所用核数的增加而增加的。
另一方面，从加速效率图中不难看出，两个算例在所用的核数$M \leq 10$时都有所下降，
但是当所用的核数$10 \leq M \leq 16$后，加速效率慢慢趋于稳定，稳定值均大于$0.5$。

%值得一提的是，比较$N=10$与$N=20$这两个算例的并行结果，我们发现当网格划分越精细时，
%即所需运算量越大时，并行的效果越好。
%这是因为，当运算量不大时，核与核之间数据交换所用的时间与核内的计算时间相比，占更大的比例。
%而当运算量增大时，核与核之间数据交换所用的时间增加的慢，核内的计算时间增加的快，结果是上述比例降低，
%也就表现为并行效率更好。
\begin{figure}[h!]
	\center{
  	\includegraphics[width=0.48\textwidth]{files/accel_ratio.jpg}
	\includegraphics[width=0.48\textwidth]{files/accel_effects.jpg}}
  	\caption{并行程序的加速比与并行效率}
	\label{accel}
\end{figure}



\section{总结}
通过这次作业，我学习到了很多很多。

（1）因为采用的是显式计算格式，所以在程序编写方面省了不少时间。
然而由于时间步长的限制，计算充分长时间后温度的分布所需要的计算时间比较长，即算的很慢。
如果采用无条件温度格式，将会大大提高计算效率。

（2）事实上，在并行策略上可以使用的方法多种多样。
然而由于我们需要解决问题的对称性与边界条件的简单性，我们采取了最最简单平均分配网格到各个核的方法。
另一方面，由于我们将程序放到cluster上运行，可以使用的核数较多，所以我们的并行加速比基本呈线性分布；
并行效率也在1左右。

（3）通过在不同平台上编译程序的经验，学习到了很多Linux系统下的知识。
比如配置PATH环境变量，安装软件并解决依赖树问题，查看CPU信息等等。

（4）通过程序运行结果的分析，我感受到了数据传递对并行效率的影响是很大的，应该尽量避免数据传递的次数
以及传递数据的大小。此外，应该尽量避免在程序中添加同步操作——强制各个进程同步地开始下一步
计算。在我们的平均分配网格到各进程的秉性策略中，这是不需要的，只会影响运行时间。



\begin{thebibliography}{99}
\bibitem{pre1} 李荣华，刘播，《微分方程数值解法(第四版)》，高等教育出版社，2009。
\bibitem{pre2} 张林波，迟学斌，莫则尧，李若，《并行计算导论》，清华大学出版社，2006。
\end{thebibliography}









\end{document}


